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C^ Abstract We study systems of close orbiting planets evolving under the influence 

^S| of tidal circularization. It is supposed that a commensurability forms through the 

j—{ action of disk induced migration and orbital circularization. After the system enters 

3 an inner cavity or the disk disperses the evolution continues under the influence 

' of tides due to the central star which induce orbital circularization. We derive 

Cn approximate analytic models that describe the evolution away from a general first 

^^ order resonance that results from tidal circularization in a two planet system and 

^.^ which can be shown to be a direct consequence of the conservation of energy 

CLj and angular momentum. We consider the situation when the system is initially 

Ft'] very close to resonance and also when the system is between resonances. We also 

perform numerical simulations which confirm these models and then apply them to 

two and four planet systems chosen to have parameters related to the GJ581 and 

HD10180 systems. We also estimate the tidal dissipation rates through effective 

O quality factors that could result in evolution to observed period ratios within the 

<4_i lifetimes of the systems. Thus the survival of, or degree of departure from, close 

i2 commensurabilities in observed systems may be indicative of the effectiveness of 

' ' tidal disipation, a feature which in turn may be related to the internal structure 

of the planets involved. 

^ Keywords Planet formation • Planetary systems • Resonances • Tidal interactions 

00 
00 

^— s Planetary systems containing hot Neptunes and hot super-Earths have been ob- 

^^ served recently. A system of this kind consists of the four planets around the 

,-H M-dwarf GJ 581 (Bonfils et al. 2005, Udry et al. 2007, Mayor et al. 2009a). The 
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projected masses of the planets are 1.9, 15.6, 5.4 and 7.1 M® and the periods 
are 3.15, 5.37, 12.93 and 66.8 days, respectively. Other such multiple systems are 
that around HD 40307 (Mayor et al. 2009b) which consists of three planets with 
projected masses of 4.2, 6.9 and 9.1 M® and periods of 4.31, 9.62 and 20.46 days, 
respectively and that around HD 10180 (Lovis et al. 2010) which consists of seven 
planets with projected masses 1.35, 13.10, 11.75, 25.1, 23.9, 21.4 and 64.4 M® and 
periods of 1.18, 5.76, 16.36, 49.74, 122.76, 601.2, and 2222 days respectively. The 
innermost member of the latter system has yet to be confirmed. 

Migration due to tidal interaction with the disk is a possible mechanism through 
which planets end up on short period orbits, as in sttu formation implies very 
massive discs (e.g., Raymond et al. 2008). Terquem & Papaloizou 2007 (see also 
Brunini & Cionco 2005) proposed a scenario for forming hot super-Earths in which 
a population of cores that formed at some distance from the central star migrated 
inwards due to interaction with the disk. These collided and merged as they went . 
This process could produce systems of planets with masses in the earth mass range, 
located inside an assumed disk inner edge, on short period orbits with mean mo- 
tions of neighbouring planets that frequently exhibited near commensurabilities. 
However, tidal circularization of the orbits induced by tidal interaction with the 
central star, together with later close scatterings and mergers tended to cause 
the system to move away from earlier established commensurabilities to an extent 
determined by the effectiveness of these processes. 

Papaloizou & Terquem (2010) considered the system around HD 40307 for 
which the pairs consisting of the innermost and middle planets and the middle 
and outermost planets are near but not very close to a pair of 2:1 resonances. In 
spite of this it was found that secular effects produced by the action of the resonant 
angles coupled with the action of tides from the central star could cause the system 
to increasingly separate from commensurability. Resonant effects can arise even 
when departures from strict commensurability are apparently large because tidal 
circularization produces small eccentricities which, for first order resonances, can 
be consistent with resonant angle libration (see Murray & Dermott 1999). 

In this paper we undertake a further study of systems of close orbiting planets 
evolving under the infiuence of tidal circularization. We present simple analytic 
models describing the evolution away from a general first order resonance for a two 
planet system under the infiuence of tidal circularization, describing the situation 
both when the commensurabilty is very close and also when the system is between 
resonances. We also perform numerical simulations of two and four planet sys- 
tems chosen to have parameters related to the GJ581 and HD10180 systems . We 
consider the situation when various commensurabilities result through the action 
of assumed disk induced migration and orbital circularization rates, estimating 
the magnitudes of tidal quality factors that could produce evolution to observed 
period ratios within the lifetimes of the systems. 

The plan of this paper is as follows. In sectiorJ2?T]we consider a coplanar system 
of planets in near circular orbits in which orbital energy is dissipated while its total 
angular momentum is conserved. The system is expected to spread in a similar 
manner to a viscous accretion disk (see Lynden-Bell 8z Pringle 1974). For a two 
planet system, this radial spreading will always lead to evolution away from an 
initially close commensurability. 

In sections |2.2[ |2.3| and [3J we carry out an anlaytic study of two planets near 
to a first order commensurability under the infiuence of tidal circularization. For 
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planets in the mass range we consider, it is readily estimated that tides raised 
on the planet are very much more important than tides raised on the star (eg. 
Goldreich & Soter 1966, Barnes et al. 2009). In addition the orbital decay timescale 
due to tides raised on the star may be estimated to be much longer than any 
timescale of interest (eg. Barnes et al. 2009) . Thus tides raised on the star have 
been neglected. 

We consider the initial evolution away from a close first order commensurability 
in section |3.1| and go on to consider the case of evolution when the commensura- 
bility is not so close, or the system is between commensurabilities in section [3. 2 1 . 
As expected from the simple arguments given in section |2.1[ the system departs 
from an initially close commensurability moving to a neighbouring one. 

We go on to perform numerical simulations of multiplanet systems in section 
|4] Various commensurabilities between pairs of planets are set up by applying 
dissipative forces assumed to arise from a disk, that lead to orbital migration and 
circularization. These forces were then removed corresponding to assumptions of 
either entry into an inner cavity or removal of the disk. The evolution of the system 
under tidal circularization caused by interaction with the central star was then 
followed. For illustrative purposes we consider two planet systems with parameters 
corresponding to the two innermost planets in the GJ581 system. In section [4.5| we 
consider a system that formed a 3:2 commensurability which then evolved under 
orbital circularization indicating that the model system could attain the period 
ratio appropriate to the actual system if the tidal parameter Q' introduced by 
Goldreich & Soter (1966) ~ 100. The situation when the system began with disk 
parameters that led to a 5:3 commensurability is then similarly studied in section 
|4.6[ We go on to consider the effect of adding the additional planets in the GJ581 
system in section [4?7l 

As there are examples of low mass planetary systems such as HD 10180 which 
have separations of pairs of planets, the third and fourth innermost in that case, 
that indicate there may have been a past proximity to a 3:1 commensurability, we 
consider an exploratory simulation of a system for which the initial disk evolution 
sets up a 3:1 commensurability in section [4.8[ Finally in section [5] we summarize 
and discuss our results. 



2 Commensurabilities and tidal circularization in planetary systems 

We begin by considering the evolution of planetary systems undergoing tidal cir- 
cularization in a general way and then move on to consider simple analytic models 
of two planet systems that can be close to first order commensurabilities. 



2.1 Two interacting planets in circular orbits for which energy is dissipated at 
fixed total angular momentum 

Consider two interacting planets with orbital energies E\ and E-2 respectively. The 
associated orbital angular momenta for assumed circular orbits are —2Ei/ni and 
—lEilrii respectively. Here n\ and 712 are the mean motions associated with the 
two planets. Suppose now the system dissipates energy while conserving its total 
angular momentum. This is expected to be the case during orbital circularization 
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when this occurs as a resuh of steUar tides dissipated in the planets because the 
planets themselves cannot contain a significant amount of angular momentum. 
Accordingly we write 

dEi dEh _ _„ ,,^ 

dt ^ dt ~ ' ^ ' 

where C is the rate of energy dissipation. Angular momentum conservation implies 
that 



from which we obtain 



and 



1 dEi 
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1 — 712 /n\ 
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1 — n\ln2 



(2) 



(3) 



(4) 



Supposing that n\ > n2, the above two equations imply that planet 1 moves 
inwards losing energy while planet 2 moves outwards, taking up the angular mo- 
mentum lost by planet 1. This is the generic form for the evolution of an accretion 
disc (see Lynden-Bell & Pringle 1974). 

We now go on to discuss some simplified models for the interaction of two 
planets that may be either very close to or some distance away from a strict first 
order commensurability. In these contexts we show how tidal dissipation induced 
by forced eccentrcities causes the system to separate. When the resonance is close, 
this causes the system to depart further from commensurability. 



2.2 Coordinate system 

We consider a general system of A*' planets orbiting a central mass. We adopt 
Jacobi coordinates (Sinclair 1975, Papaloizou & Szuszkiewicz 2005) for which the 
radius vector of planet i, Tj, is measured relative to the centre of mass of the system 
comprised of a dominant central mass M and all other planets interior to i, for 
i = 1,2..., N. The planets are assumed to maintain an ordering with increasing 
i corresponding to greater distances \ri\ from the dominant central mass. Thus 
the innermost planet has i = 1. The Hamiltonian, correct to second order in the 
planetary masses, can be written in the form: 

N , 

2 GM,m, 






Here Mi = M -\- rrii and vij = r^ — r, 
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The equations of motion for motion for planet i assumed to move in a fixed 
plane, about a dominant central mass, may be written in the form (see, e.g., 
Papaloizou 2003, Papaloizou & Szuszkiewicz 2005): 



(6) 

(7) 

(8) 

(9) 

Here the orbital angular momentum of planet i which has reduced mass m^ = 
ithqM / {M + mio), with rriiQ being the actual mass, is Li and the orbital energy is 
Ei. For motion around a central point mass M we have: 
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, = m,^GM,a,{\-el), (10) 

S, = -^, (11) 

where Mi = M + rriiQ, ai denotes the semi-major axis and e^ the eccentricity of 

planet i. 

The mean longitude of planet i is Aj = nj(t — t^i) -\- vji, where rn = ^GMi/a^ 

is its mean motion, with foi denoting its time of periastron passage and vji the 

longitude of periastron. 

From equations ([6| and ([7| an equation for the evolution of the eccentricity of 

planet i may be readily obtained in the form 

- = ^^fl?f\/W-l)+|^). (12) 

eimiUiaf \d\i \\ J du^i J 

The Hamiltonian may be expanded in a Fourier series involving linear combina- 
tions of the (2A'' — 1) angular differences -Wi — zuiji = 2,3..N and Xi~-cui,i = 1,2, ..A''. 
In the limit of small eccentricities of interest here, only terms that are of first order 
in the eccentricities need to be retained (terms that are of zero order do not lead 
to changes to eccentricities or to resonances) . If this is done the possibility of first 
order resonances, for which the ratio of the periods of two planets is the ratio 
of successive integers, is allowed for. The above approximation scheme should be 
valid when circularization times are small enough to ensure that the eccentrici- 
ties remain small. This situation is realized for examples of low mass protoplanets 
migrating in protoplanetary discs (Papaloizou & Szuszkiewicz 2005). 

Near a first order p -|- 1 : p resonance, p being an integer, we expect that terms 
in the Hamiltonian involving angles of the type 4'p,j,i,j = (p+ l)Aj — pAi — vuj, and 
4'p,j,i,i = (p+ l)Aj — pAi ~rx7i, where the subscripts on the left hand side correspond 
to those on the right hand reading from left to right, will be slowly varying and 
thus be dominant. Accordingly we shall retain only terms of this type. Motion 
away from resonances may also be considered having made this approximation 
although neglected high frequency modulations may be more significant then. 



The Haniiltonian may be written in the form 

N N-1 N 



^ = E^^ + E E^^^ 
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(13) 



i=l j-i+1 



where the component of the interaction Hamiltonian Hij that is first order in the 
eccentricities is given, given that j > i, by 



Gmim 



Hi,j = > {^jCp,j,i,j cos(0pj,ij) + eiDpj,i^i cos(0pj,i,i)) , (14) 



with 



dih^^Ux)) 



^^■^•^■^' ~ 2 r dx 



+ (2p + V)h^^1^ (x) - ^x&l and 



(15) 



D„ 



l(.^!&.....,-(. 



(16) 



Here tif} Ax) denotes the usual Laplace coefficient (e.g. Brouwer & Clemence 1961) 

with the argument x = ai/aj and Sf denotes the Kronnecker delta. We remark that 
the subscripts associated with the coefficients Cp^^ij and Dp^^i^ correspond to the 
related angles as in ( 14 1 . We shall also make the approximation of replacing Mi 



by M and equivalently ttijq by m^. 

The governing equations for motion, retaining only terms that are of the low- 
est order in the eccentricities, follow from Hamilton's equations ^- ^ for the 
Hamiltonian ( 13 1 discussed above as 
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~dt 
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(17) 



oLrTrijp I eiJJpji_iSui[(ppji^i) ejUp^j_ijSin.y<ppj^ij) 



dzui 
"dT 
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E 
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36*771^(^-1-1) f ejDp,i^jjsm{(j)p^ij^j) eiCp^ij^iSin{<j)p^i^j,i) 



+ 



EGmjDpj^i^iCos{<j)pj.i^i) ■^-^ GmjCp,ij,iCos[(j>p^ij^i) 



j=-j+l 



i=i 



(18) 
(19) 



In addition, consistent with the above approximation scheme, the rate of change 
of the mean longitudes may be obtained from 



dXi 

—— = riA 
dt ' 

which also enables evaluation of the rate of change of the angles . 



) J,2,Z! T^P,J^i,3 



(20) 
etc. 
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2.3 The incorporation of disk tides 

We incorporate the effects of orbital circularization by adding additional terms to 



the right hand sides of equations ( 17 1 and ( 18 1 . Equation ( 17 ) is modified through 
the straightforward prescription 



dt 



dcj 
dt 



(21) 



where t^^t is the circularization time for planet i. Similarly equation ( 18 ) is modified 
according to 

+ ^. (22) 



dui 
~dt 



drii 
~dt 



ta 



This adjustment is necessary to account for the orbital energy dissipation occurring 
as a result of circularization correct to the lowest order in e^. This dissipation is 
assumed to occur with out changing the angular momentum of the system because 
the planets can only potentially contain a negligible amount of angular momentum 
compared to that in the orbit. It follows from the energy dissipation rate for planet, 
i, given by 



dt 



{^-e1)tc,i 



(23) 



For small eccentricities, e^ may be neglected in the denominator of the above 
expression and the total rate of energy dissipation in the system is obtained by 
summing over all planets. 



3 Two planets in a p + 1 : p commensurability 



It is possible to investigate solutions of equations (171 - (191 modified to incor- 



porate circularization, that illustrate the geometrical separation of the system as 
energy is dissipated while the total angular momentum is conserved, in a number 
circumstances. 



3.1 A tight commensurability 



We begin with an example where two successive planets k and fc + 1 maintain 
a p + 1 : p commensurability with the associated angles in a state of at most 
small amplitude libration while their semi-major axes separate. We later go on 
to consider a simple restricted example where the angle circulates. The effects 



of planets other than the resonant pair is neglected. Equations( 17l - (191 with 



the modifications given by (21 1 and (22) to incorporate circularization give the 



governing equations for planet k in the form 
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dt al \ flfc+i 

ek+lCp^k + l,k,k + lS^^{(t>p,k + l,k,k + l)\ Srtfeej, 



O-k + l / tc^f^ 

dwk _ Gmj,j^-i^Dpi,j^i^j,j,cos{<j)p^k+i,k,k) 
dt fifc'^fc'ifc'Jfe+i 

Similarly the governing equations for planet fc + 1 are given by 



rfwfc+i _ 3G'mfc( p+ 1) /^e^Dpfc+i fc fesin(0p fc_|_i fe fc) 



(26) 



(27) 



e*;+iC'p.fe-|-i.fe.fe-|-i sin(0p.fe_|_ijejt_|_i) \ 3nj,+ie^,^j^ 

H (,28j 



(29) 



i^fe+l / '^c.fe + l 

du^k+i _ Gmi,Cp^k+i,k,k+iCos{(l)pi,^ij,k_^.i 
dt ek+iUk+ial^^ 



Setting (t>p^k+i,k,k -^ 'l>p,k+i,k,k ± ^'t>p,k+i,k,k and (pp^k^i^k^k^i -> 0p^fe+i^fc^fe+i ± 
^</'p,fc+i.fe,fc+i) where the positive sign is taken when the equilibrium value of the 
angle is zero and the negative sign is taken when it is tt, and A indicates a small shift 
such that the sines of the angles may be replaced by the angles themselves. Then 
assuming that the evolutionary time scale is much longer than the circularization 
times so that the time derivatives of the eccentricities may be neglected, we can 
then find expressions for the small angular shifts in the form 

. , _ ekUkalak+i .„„^ 

A(pp,k+iM,k - --7^ ^ 7 — ■ (■iU) 

., ek+ink+ial+i 
A(pp^k+i,k,k+i - --PT—-F, ; • (■^J-j 

Substituting these into the equations for the evolution of the mean motions yields 

ririfc ^ 3(p+l)egnfc 3pmfc+ieg^;^Wfc_nag_^^ 
dt tcM mkultc^k+i 



"k 
dnk+i ^ 3(p+ l)mkelnkal _ 3pe|_)_^nfc_|_i 



(33) 



The above pair of equations express the conservation of energy and angular mo- 
mentum for the system in the limit of small eccentricity. We remark that the latter 
follows in the form 

2 duk^i 2dnk , , 

rrik+iak+i^j-^ 1- "ifca^.-^ = 0, (34) 
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while the former follows from using the fact that Ej oc rnjU-' to find equations 
for dEj/dt,j = k,k + 1 and then adding. We may also obtain an equation showing 
how the period ratio increases with time in the form 



d 
dt 



^k+i 



Sn^J 



+ 



p4+i 



ic,kJk+l ic,k+lJk 



(35) 



where J^ = mj^aj^rn^, and J = J}^-\- Jk+\- In order to proceed further we need to 
calculate the eccentricities. These may be obtained from the governing equations 



for the evolution of the angles that may be obtained from (201, (261 and ( 29 1 in 
the form 



0'Vp,k+l,k,k 
Jt 

t>p,k+l,kM+\ 

dt 



(P + l)j^fc + l - VTT-k - 



[P + l)j^fc+i - vnk - 



Gmk+iDp^k+l,k,k COs{(l)p^k+l,k,k) 

GmkCp^k+i,k,k+i<^os{(l>pk^ikf^^i) 
efc+i'^/c+i'ifc+i 



(36) 



(37) 



As the angles are quasi-steady and close to zero or tv, these expressions enable 
the calculation of the squares of the eccentricities ej, and ej._|_i which are required 



in order to calculate the rate of period separation through ( 35 I . They are found 
to be given by 



2 



2 
Cfc+l 



Gmk+iDp^i^^i^k,k 



nkO-lak+iiiP + l)wfe-|-i - puk] 

GmkCp^k+i,k,k+i 
■n-k+ial^i [{P + l)n-fe+i - P"-/c] 



and 



respectively. 



(38) 
(39) 



Using these in (351 we obtain 



A. 
dt 



nk 



P+1 



n-k+i 



(40) 



where 



(p+l) /Gmj._,_i_Dp fc_|_i fe fe 



tc,kJk+i \ prikTik+iaj^ak+i 



ic,k + lJk 






(41) 



When the system starts to move away from a commensurability taken to be exact 
at i = 0, we may treat the right hand side of (401 as being constant and integrate 
with respect to time to obtain 



^k 



P+1 _ hnkJ 
P \nk+i 



Ft 



1/3 



(42) 



A similar scaling for which the separation from a commensurability increases oc t^'^ 
was obtained for a three planet system by Papaloizou & Terquem (2010). 
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3.2 The interaction between two planets away from a close commensurability 



In this case we again assume interaction between planets k and k + 1. In this case 
we consider the situation away from a strict commensurability where significant 
libration or circulation may occur. This is a natural development as tidal evolu- 
tion causes the system to evolve away from a tight commensurability of the type 
described above towards such a situation. We make the additional simplification 
of assuming that rrifc^i ^ ruf;. In that case, a circular restricted 3 body prob- 
lem may be adopted. Only the motion of planet k is considered with ej.^i = 0. 
Equations (17 1 - (201 apply and as ej^_^_i = 0, only terms involving the angles 
4'p,k+i,k,kjP = 1,2... appear. These give the equations governing the evolution as 



rffifc _ 

dt 


fifc 

tc,k 


duk 
dt 


3nkel 

tc,k 


d<j)r,k+l,k,k 


(r + 1) 


dt 


r = 1,2,3 . 





E 



Gmf;^iDp^k+i,k,k s>in{(l>p^k+i,k,k) 



■\r^ SGmfc+ipefcDp^fe+i^fc^fe sin(<^p^fc+i^fc^fc) 



p=i 



afc+iflfc 






ek^ka/^ak+i 



(43) 



(44) 



(45) 



Although we consider the effect of more than one angle, we focus on a particular 
one with r = q which might be considered to be the one closest to resonance, though 
that is not essential. Setting x = ef;Cos{(j)q]^^ikk) ^nd y = ekS^'^{4'q,k+i,k,k) in 
equations (431 and ( 45 ) leads to a system that, unlike the original one, does not 
contain an apparent singularity as ej, — > in the form 



dx 

'dt 



'^qV- 



-T 2^ a?pSin[(p-g)(Afe+i - Afe)] 



ic.k 



p=ljtq 



^ = -a 

dt 11 ^^ 



(46) 



+ u}qx--^- V apCos[(p-g)(Afe+i - Afc)], (47) 

p=l^q 



where, recalling that k is fixed, we define aq = Gmfe^iDq j^^i fcjj,/(rn;a^aj,^i) 
and ujq = {q + l)nk+i - qn^. 

We now remark that if we consider the limit ej, — >■ 0, which occurs far enough 
away from resonance, we may neglect the evolution of nj. and asume that it remains 



constant and equal to nj-g. To see this it follows from (431 that e^. scales as mfc^i 
while (44 1 then indicates that the change in n^., 5rn. = nj. — n^-g scales as mj.,^ or 



e^.. Accordingly, for q of order unity, the induced variation of Wq, Suiq, is such that 
Sujq ~ nfcoe|. By comparing the variation of the first two terms on the right hand 
side of equation (47), it readily follows in the low eccentricity limit that provided 
'T'feoefc < Og or 



f Gm^.^iDqj.^11; J. 

fife < 2 2 



1/3 



(48) 
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the variation of nj. may be neglected so that it may be taken to be equal to 
Uf^Q. Similarly aj. is replaced by the corresponding fixed value a^-o- The above 
approximation scheme applies in the low eccentricity limit or sufficiently far away 
from strict commensurability such that \(jjq\ = \{q + l)ni.+i — qrikol ^ ^kn '^i • 



Given that in the same approximation ( 20 1 implies that 



Xk+i — Afe = {rik+i — rn.o)t, equations (461 and (47) describe a linear system with 



prescribed harmonic forcing that is easily solved exactly. The solution in the limit 
tc ^- oo may be written 



. = .0+ y "p'^«^[(P-^)(^^+i-^fc)l (49) 

■^ (P + ^)nk+i - pnk 

_ y^ QpSin[(p-g)(Afc_n-Afc)] 
^ {p + l)nk+i - puk 

where xq = Uq/uug. This indicates oscillation about the mean value of a; = xq. 
In the absence of the periodic forcing circularization would cause the solution to 
approach x = xo,y = corresponding to a precise commensurability with zero 
libration amplitude. When the forcing is present, there is either libration or circu- 
lation depending on the ratio of the forcing amplitude to xq. When g + 1 : q is the 
closest commensurability, this can be small resulting in small amplitude libration. 
When some other commensurability is dominant, the motion in the {x,y) plane is 
around an approximately circular curve that encloses the origin and so corresponds 
to circulation. Thus as the system moves through a commensurabilty the motion 
is expected to change from circulation to libration to circulation (see Terquem & 
Papaloizou 2007 for an example of evolution away from a first order commensu- 
rability driven by orbital circularization) . Note that in all of these cases the time 
dependent averages of quantities such as Ck cos{(l>q,k+i,k,k) = ^ ^^nd cos{(l>q,k+i,k,k) 
are generally non zero (see also Papaloizou & Terquem 2010). 

In order to calculate the rate of energy dissipation resulting from orbital cir- 
cularization and hence the rate of evolution it causes, we require the time average 
of the square of the eccentricity. This is given by 

oo 2 

(el) = {^' + y')=T. u ^u "' P- (51) 

^ [{P+ l)'^fc+l -pnkr 
We remark that this expression connects to that found for the tight resonance 



example (381. The latter expression is identical to that given by (51) if only one 
term is retained that corresponds to the tight commensurability considered. Thus 
we expect the evolution to continue to be dominated by the closest comensurability 
until another becomes closer and takes over governing the evolution. 



3.3 The orbital evolution of the planet 

The rate of change of the orbital energy may be obtained from consideration of 
(51 1 and (23) together with the discussion leading to equation (pi given in section 
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I2.1l with the result that 

2 / \ 2 2 ^^'^ 2 



dt 3 dt\nk+ij {l-nk+i/nk)tcM "^^[{p+^hk+i-pn-kV 

(52) 
which means that the orbit of mj. contracts and separates from that of 7nj._|_i. In the 
limit m,/;/mj._|_i — >■ in which n^j^i becomes fixed, and close to a commensurability 



( 52 I becomes equivalent to ( 40 1 . Thus it enables the discussion of the situation 
corresponding to a tight commensurability to be extended to conditions away from 
close commensurability. In view if the fact that n/j/nj. _|_i must always increase, 
such a discussion leads to the conclusion that the evolution will be controlled by 
successive first order comensurabilities as the system widens (see also Terquem & 
Papaloizou 2007). 



4 Numerical Simulations 

We here describe simulations of model planetary systems in which commensurabil- 
ties have been formed subsequently evolving under the infiuence of circularization 
tides. 



4.1 Model and initial conditions 

We consider a primary star together with N planets embedded in a gaseous disk 
surrounding it. The planets undergo gravitational interaction with each other and 
the star and are acted on by tidal forces from the disk and star. The system is 
solved as an A^-body problem. Tidal interactions are incorporated by applying 
appropriate dissipative forces (see Terquem & Papaloizou 2007 and Papaloizou & 
Terquem 2010 for more details and examples). The equations of motion may be 
written as 



9 AT „ , 

d Fj _ GMr^ Y^ Gmjo (r^ 



di2 |r,|3 ^ ir,-r,|3 

3 = 1^1 



-r + r, + rr , (53) 



where Af , rrijQ and tj denote the mass of the central star, that of planet j and the 
position vector of planet j, respectively. The acceleration of the coordinate system 
based on the central star (indirect term) is given by 

T^. (54) 

and that due to tidal interaction with the disk and/or the star is dealt with through 
the addition of dissipative forces (see Papaloizou & Larwood 2000). Thus 

r 1 *i 2 fdr, \ 2 f dvi \ . 

tmg,i dt Iril-^te^i \ dt J ti^i \ dt J 

where tmg,i, te.i and ti^i are the timescales over which, respectively, the angular mo- 
mentum, the eccentricity and the inclination with respect to the unit normal e^ to 
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the assumed fixed gas disk midplane change. Evolution of the angular momentum 
and inclination is assumed to be due to tidal interaction with the disk, whereas 
evolution of the eccentricity is assumed to occur due to both tidal interaction with 
the disk and the star. We have: 

— = ^^ H . (56) 

f . j-d f . ' ^ ' 

'* c i '* 

where t^ ^ and ic,i &re the contribution from the disk and tides raised by the 
star, respectively. Relativistic effects are modeled through Tr ( see Papaloizou & 
Terquem 2001). 

Because a low mass planet cannot contain a significant quantity of angular 
momentum, tides raised on it by interaction with the star are assumed not to 
modify the angular momentum of the orbit. We remark that the orbital decay 
timescale, due to tides raised on the star, is readily estimated to be much longer 
than any timescale of interest (eg. Barnes et al 2009) thus these tides are ignored 
from now on. 



4.2 Orbital circularization due to tides from the central star 

The circularization timescale due to tidal interaction with the star, in the small 
eccentricity limit appropriate here, is taken to be ( Goldreich & Soter 1966) 

'— «'(^)"(S)"(S)"«'>-- *-) 

where a^ is the semi— major axis of planet i. Here we have adopted a mass density of 
1 g cm for the planets (uncertainties in this quantity could be incorporated into a 
redefinition of Q'). The parameter Q' = 3(5/(2fe2), where Q is the tidal dissipation 
function and ^2 is the Love number. For solar system planets in the terrestrial 
mass range, Goldreich & Soter (1966) give estimates for Q in the range 10-500 
and fc2 ~ 0.3, which correspond to Q' in the range 50-2500. We remark that this 
parameter should be regarded as being very uncertain for extrasolar planets. As 
computations with increasing Q' become prohibitive on account of long evolution 
times, we have considered values of Q' of 1.5 and 3 in this paper. However, we 
have obtained scaling relations which indicate how to scale results to larger Q' . 



4.3 Type I migration 

When a planet is in contact with the disk, disk-planet interactions occur leading to 
orbital migration as well as eccentricity and inclination damping (e.g.. Ward 1997). 
However, the migration rates to be used are uncertain even when the disk surface 
density is known, largely because of uncertainties regarding the effectiveness of 
coorbital torques (e.g., Paardekooper & Melema 2006, Pardekooper & Papaloizou 
2008, 2009). In this context there are indications from modelling the observational 
data that the adopted type I migration rate should be significantly below that 
predicted by the linear calculations of Tanaka et al. (2002) (see Schlaufman et al 
2009). Hence we have carried out simulations with tmg.i and tf j for any system 
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Fig. 1 The evolution of two planets that form a 3:2 commensurability. The early evolution 
of the period ratio during convergent migration is shown in the upper left panel. The upper 
right panel shows the evolution of the period ratio under orbital circularization after disk 
migration ceases. The uppermost curve is for Q' = 1.5 and the lower curve is for Q' = 3. The 
triangles/diamonds correspond to the analytic predictions made from equation \42\ adapted 
to the case of a 3:2 commensurability for Q' = 1.5/Q' = 3 respectively. The evolution of 
the eccentricity of the outermost planet is plotted in the lower left panel for Q' = 1.5. The 
evolution of the resonant angle 3A2 — 2Ai — -uui is plotted in the lower right panel for Q' = 1.5. 



taken, as for type I migration, to be proportional to l/nij and adopted f j j = 
if j. A range of scaling constants was explored. These are quoted together with 
corresponding numerical results below. 

We remark that provided that eccentricity damping limits eccentricities to 
small values, the commensurabilities that are formed in the system as a conse- 
quence of convergent migration depend on the ratio of the adopted migration rate 
to the local orbital frequency, with commensurabilities of low order and low degree 
forming when this ratio is small. 
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4.4 Numerical results 

4.5 A system with a 3:2 commensurability 

For the calculatioas presented in this section we adopted masses for the two planets 
and the central star that coincided with those for the star and two innermost 
planets of the GJ581 system. Thus the inner planet was taken to have a mass 
m\ = 1.94Af0 and to be in circular orbit at Q.\&au. The outer planet was taken 
to have a mass rrii = 15.64Mop;„c, and to be in a circular orbit at 0.32au. Tests 
indicate that the results of simulations of the type described here do not depend 
on the longitudes at which the planets are inserted on such circular orbits. The 
central mass was O.SIM©. The initial semi-major axes were chosen to be larger 
than the corresponding ones in the GJ581 system so as to allow for some inward 
migration. The disk migration and circularization rates adopted were given by 

t„,„ = 4.375 X lO'^— w. and t^ » = 5 x 10^— yr. (58) 

However, they were only applied when the semi-major axis of a planet exceeded 
0.041aM. This procedure results in the final semi-major axis of the outer planet 
to coincide with the second planet in the GJ581 system. The termination of disk 
migration could be regarded as either being due to entry into an inner cavity, or 
simply removal of the disk. The migration rate was chosen so as to enable the 
planets to settle into a 3:2 commensurability through convergent migration. A 
very much slower rate would allow trapping in a 2:1 commensurability, while a 
very much faster one would result in the system passing through the 3:2 commen- 
surability (see eg. Papaloizou & Szuskewicz 2010). We remark that although the 
specific parameters chosen correspond to the GJ581 system, the aruments pre- 
sented above indicate that the form of evolution we find should be generic for 
two low mass planets attaining a first order commensurability through convergent 
migration. 

The evolution of the system is illustrated in Fig. [l] The early evolution of the 
period ratio during convergent migration is shown in the upper left panel. It is seen 
that the system is trapped in a 2:1 commensurability for a while before escaping 
to be subsequently trapped in a 3:2 commensurability. After about 2 x W^yr. the 
forces from the disk cease to act and the system evolves under tidal circularization. 
The upper right panel of Fig. [l] shows the evolution of the period ratio. Results 
for simulations with Q' = 1.5 and Q' = 3 are illustrated and compared to analytic 
predictions derived from equation ( |42[ ) adapted to the cases on hand. Interestingly 
the numerical results are in quite good agreement with what is expected from the 
analytic discussion given in section |3] which assumed a small libration amplitude 



and which led to equation ( 42 I , even in regimes where the amplitude of libration of 
the resonant angle is quite large. However, the simulations show additional sudden 
small jumps in the period ratio which occur when the system passes through the 
5:3 resonance. This jump was larger for the Q' = 1.5 case than for the Q' = 3 case. 
The evolution of the resonant angle 3A2 — 2Ai — -cui for Q' = 1.5. shown in Fig. [l] 
indicates an increasing amplitude libration that tends to break down near the end 
of the simulation when the period ratio ~ 1.7 as in GJ581. But note that there 
is also a short temporary breakdown as the system passes through 5:3 resonance. 
Note that the anlaytic treatment suggests the time for the period ratio to evolve 
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Fig. 2 The evolution of two planets that form a 5:3 commensurability. The early evolution of 
the period ratio during convergent migration is shown in the upper left panel. The upper right 
panel shows the evolution of the period ratio under orbital circularization after disk migration 
ceases. The uppermost curve is for Q' = 1.5 and the lower curve is for Q' = 3. The evolution 
of the eccentricity of the outermost planet is plotted in the lower left panel. The evolution of 
the resonant angle 3A2 — 2Ai — vo\ is plotted in the lower right panel. 



from 1.5 to 1.7 to be ~ 5 x lO^j/r. The simulation with Q' = 1.5. rather fortuitously 
agrees very well with this. The analytic prediction for Q' = 3 is lO^j/r. while the 
simulations discussed in this and the next section indicate 1.3 x lO^yr. Given the 
expectation that evolution times are oc Q' , this indicates that values of Q' as large 
as a few hundred could have allowed the period ratio to move from 1.5 to the 
present value within the lifetime of the system. 



4.6 A system with a 5:3 commensurability 



For the calculations presented in this section we adopted the same values for the 
central mass and the planet masses as in section [ITsl However, we adopted initial 
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conditions, migration and circularization rates so as to enable the system to settle 
into a 5:3 commensurability. Thus the inner planet was started in circular orbit 
at O.OSau and the outer planet in a circular orbit at 0.16aw in this case. The disk 
migration and circularization rates adopted were given by 

t^^ = 1.75 X W^^yr. and t^ = 2 x W^^yr. (59) 

mi ' TUi 

Thus the convergent migration rate was two and a half times faster and the circu- 
larization rate four times slower than for the calculation in section |4.5| However, 
they were applied in the same way. The faster migration rate and the slower ec- 
centricity damping rate allows trapping in the 5:3 resonance. 

The early evolution of the period ratio during convergent migration is shown 
in Fig. [2J It is seen that the system becomes trapped in a 2:1 commensurability 
before escaping to form a 5:3 commensurability. After about 8 x lO^yr. the forces 
from the disk cease to act and the system evolves under tidal circularization. 

The upper right panel of Fig.[2]shows the evolution of the period ratio under or- 
bital circularization after disk migration ceases for Q' = 1.5 and Q' = 3. Although 
the system started in a 5:3 commensurability, the evolution can be regarded as 
matching onto that illustrated in the previous section which can be regarded as 
being driven by the 3 : 2 comensurability. This is also confirmed by the evolution 
the resonant angle 3A2 — 2Ai — tui also plotted in Fig[2] We also remark that the 
time for the period ratio to move from 5/3 to 1.7 is about twice as large for Q' = 3 
as for Q' = 1.5 as expected. However, these times are only approximately 1.8 x 10® 
and 3.6 x 10® yr. respectvely indicating that values of Q' up to lO'^ could be effective 
within the lifetime of the system. 



4.7 Adding additional planets 

Here the effect of adding additional planets to the simulation described above is 
investigated. To do this we take the calculation of section [4!6| at the point at which 
forces arising from the disk cease to act. Two additional planets of masses 5.36M0 
and 7.O9M0 are added in circular orbits with semi-major axes O.QIau and 0.22aM 
respectively. These correspond to the two outermost planets in the GJ581 system. 
We remark that the eccentricities of these planets were determined to be consistent 
with zero by Vogt et al. (2010). As before we considered runs for which Q' = 1.5 
and Q' = 3. In this case the same value of Q' was adopted for each planet. 

The results are plotted in Fig. [3l The evolution in this case is for the most part 
similar to that illustrated in Fig. [2] for two planets. In particular approximately 
the same time is taken for the period ratio for the innermost pair of planets to 
move from 5/3 to 1.7. However, a significant difference is that the evolution of the 
period ratio slows down briefiy between 1.5 x 10^ and 2.0 x lO^yr. in the simulation 
with Q' = 3. During this time the eccentricity of the second innermost planet 
is increased. Although the reasons for this are unclear, it is associated with an 
interaction between the second and third innermost planets. The innermost planet 
continues to move inwards but the angular momentum ends up being transferred to 
the third rather than the second innermost planet. There does not seem to be any 
clear resonance associated with this. However, we comment that in a many planet 
system like this, we could consider a tension between possible interacting pairs. The 
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Fig. 3 The evolution of two planets illustrated in Fig.lslbut with two additional outer planets 
added after the disk migration phase as indicated in the text. The upper right panel shows 
the evolution of the period ratio of the innermost two planets under orbital circularization. 
The uppermost curve is for Q' = 1.5 and the lower curve is for Q' = 3. The evolution of the 
eccentricity of the second innermost planet is plotted in the upper left panel. The evolution 
of the resonant angle 3A2 — 2Ai — wi is plotted in the lower left panel. The behaviour of the 
angle between the apsidal lines of the orbits of the second innermost and innermost planets is 
illustrated in the lower right panel. 



second and third planets would separate on account of tidal circularization if the 
innermost planet were absent. Similarly the innermost pair can couple as in section 
|4.6[ In some circumstances, dependent on their masses, orbital parameters, and 
values of Q' etc., different interacting pairs may have varying levels of importance 
in the simulation. This requires a more detailed study than we have been able to 
perform at this preliminary stage that will be the subject of future work. 
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Fig. 4 The evolution of two planets in a 3:1 eommensurability is illustrated. The left upper- 
most panel shows the evolution of the semi-major axes of the two planets. The initial period of 
disk migration is short <'~ 4.5 X lO^yr. The subsequent evolution is driven by tidal circulariza- 
tion with Q' = 100 and the gravitational interaction between the planets. The upper left panel 
shows the evolution of the period ratio. This remains 3:1 for some time after disk migration 
has ceased before finally increasing as the planets separate. The lower left panel shows the 
evolution of the eccentricities of the two planets, the uppermost curve corresponding to the 
inner planet. The lower right panel shows the angle between the apsidal lines of the orbits of 
the outer and inner planets, which ultimately remains close to n. 



4.8 A system with a 3:1 eommensurability 



Finally we describe a situation in which a 3 : 1 eommensurability could be formed 
under convergent migration and then subsequently maintained. The parameters 
of this simulation were chosen to lead to a separation of pairs similar to the third 
and fourth innermost planets in the HD 10180 system for which there may have 
been a past proximity to a 3:1 eommensurability. 

In this case the central mass was taken to be IM®. The inner planet mass was 
taken to be 11.73M0 and the outer planet taken to be 25.07M©. Their initial semi- 
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major axes were 0.387aM and l.2au respecively. The outer planet was started in 
circular orbit while the inner planer was started at apocentre with an eccentricity 
e = 0.24. The disk migration and circularization rates adopted were given by 

tra^g =tc,^ = lAx 10'' ^yr. (60) 

These were applied only when the planets semi-major axes exceeded 0.29au. Note 
that this migration rate is very much lower than the previous cases so as to en- 
able trapping in the 3 : 1 resonance. The ecentricity damping rate is taken to be 
equal to the migration rate so that the eccentricities do not damp too quickly so 
enabling the 3:1 resonance to persist. Rates like these are not readily produced in 
calculations of disk planet interactions for which the planets are fully embedded. 
They may be possible if the planets are located within a wide cavity. However, this 
aspect remains to be investigated. Here we simply adopt these rates and explore 
their consequences. Because of the larger planetary masses and larger orbital ec- 
centricities in this run, it was possible to consider larger values of Q' . We adopted 
Q' = 100. 

The evolution of the two planets in this simulation is illustrated in Fig. [4j The 
planets undergo convergent migration and attain a 3:1 resonance. The eccentricity 
of the inner planet grows up to ~ 0.56. The growth ceases after ~ 4 x W^yr. 
when effects arising from the disk cease to act. After this time the planets evolve 
under tidal circularization. For about 50 million years the commensurability is 
maintained while the eccentricty of the inner planet decreases and that of the 
outer one increases. In this process angular momentum is transferred to the inner 
planet. However, this form of evolution cannot be maintained and it reverts to 
the situation where the planets separate in semi-major axis as descibed in section 
|2.1| while the eccentricities decrease. The period ratio secularly increases while the 
angle between the apsidal lines of the orbits of the outer and inner planets remains 
close to TT. 

Interestingly at least five resonance passages were seen during this later evolu- 
tionary stage. As the period ratio increased, these corresponded to the 19:6, 16:5, 
13:4, 10:3 and 7:2 resonances. They are manifested as local blips in the eccen- 
tricity evolution of both planets as shown in the lower left panel of Fig. |4] The 
resonance passages are of decreasing order with increasing time and so the conse- 
quent changes induced in the planetary eccentricities increase in magnitude. The 
fact high order resonances such as 19 : 6 were manifest in this run is because of 
the relatively high eccentricities, in particular of ~ 0.3 for the inner planet. 



5 Discussion 

In this paper we have studied systems of close orbiting planets evolving under the 
influence of tidal circularization. We considered the situation where the system 
evolved under the influence of disk tides to form a commensurability. After the 
disk tides ceased to operate, either because of entry into an inner cavity, or because 
of loss of the disk, the operation of tidal circularization caused increasing departure 
from any close commensurability as time progressed. 

In sectiorJ27T] we pointed out that a system of planets in near circular orbits is 
expected to separate on average as energy is dissipated while angular momentum 
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is conserved. This is also expected in the very similar situation of an accretion 
disk evolving under a viscosity (see Lynden-Bell & Pringle 1974). In the simplest 
case of two planets, this inevitable increasing physical separation has to lead to 
the increasing departure from any initial commensurability. 

In sections |2.2| |2.3| and [3J we developed a formalism that could be adapted 
study the evolution of two planets near to a first order commensurability under 
the infiuence of tidal circularization. This was then applied to a system with a tight 
commensurability in section [3.1[ An expression for the depart ure from commensu- 
rability, indicating this to be oc t^'^, was given (see equation (42)). The discussion 



was then extended to the situation when the two planets were not necessarily in a 
close commensurability in section |3.2| The orbital evolution of the planet in that 
case, leading to a neighbouring commensurability, was then considered in section 

EH 



In order to confirm the analytic modeling, numerical simulations were under 
taken in section |4] We were able to set up systems of low mass planets in vary- 
ing commensurabilities, depending on the strengths of the disk tides leading to 
orbital migration and circularization, with weaker tides in general leading to more 
widely separated commensurabilities. We focused on a two planet system which 
had the same parameters as the innermost two planets as the GJ581 system in 
section [4?5l This formed a 3:2 commensurability which then evolved under orbital 
circularization. This model system attained the period ratio of the actual system 
after ~ 10® yr. when Q' -^ 1. Simple extrapolation thus indicates that tidal evo- 
lution could have moved the system to the present period ratio of 1.7 from a 3:2 
commensurability if Q' ~ 100. Similarly the situation when the system initially 
attained a 5:3 commensurability was studied in section [Xb] In this case the evolu- 
tion quickly adapted to evolve as for the case with the initial 3:2 commensurability 
when that had reached the same period ratio. However, a larger value Q' would 
suffice to cause the period ratio to move from 5:3 to the observed one within a 
given life time. The effect of adding the additional planets in the GJ581 system 
was considered in section [4?7| In that case over the long term the extra planets did 
not greatly affect the evolution. However, for a brief period the third planet moved 
outwards taking up the angular momentum of the innermost planet rather than 
the second, with the consequence that the period separation rate for the innermost 
pair was slowed. Thus a pair of planets may not always evolve independently of 
others in the system, a feature that requires further study. 

Finally the evolution of a system that formed a 3:1 commensurability was con- 
sidered in section |4.8| The model system adopted had similar parameters to the 
third and fourth innermost planets in the HD 10180 system. This case required 
slow disk migration and weak circularization with the result that the resonance 
involved high eccentricities. Because of these the commensurability could be main- 
tained for a while under orbital circularization. However, eventually the system 
increasingly departed from it as in the other cases. Finally all of our results indi- 
cate that if Q' <~ 100, commensurabilities would have been significantly affected 
by tidal effects related to orbital circularization. Thus the survival of close com- 
mensurabilities in observed systems may be indicative of the presence of large Q' 
values, a feature which in turn may be related to the internal structure of the 
planets involved. 
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